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COMPUTER SIMULATION OF SURFACE 
AND FILM PROCESSES 


This program has three main parts: 

I. Film Formation and the Morphology of Solid Surfaces and Overlayers. 

II. Diffusion of Atoms and Microclusters at Substrate Surfaces. 

III. Crack Formation and Propagation Phenomena. 

All three investigations are based on computer simulation techniques and 
the objective is to analyze various processes at atomic levels. The progress 
made during the last six-month period is presented in this report in three 
different self-contained chapters. 
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I. FILM FORMATION AND THE 


MORPHOLOGY OF SOLID SURFACES AND OVERLAYERS 
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Objective 

The objective of this study is the analysis of surfaces and interface 
properties at microscopic levels. Two types of simulations have been 
initiated, (a) parametric investigations and (b) a specific case study. The 
parametric investigations are made for the purpose of analyzing the general 
behavior of different systems, comparatively, and to classify the important 
factors which influence both the structure and the energy-related properties 
of surfaces and interfaces for these systems. 

In this study, a molecular dynamics technique based upon Lennard-Jones 
type pair interactions is used to investigate time-dependent as well as 
equilibrium properties. The specific case study deals with systems containing 
Si and 0 atoms. In this case a more involved potential energy function (PEF) 
is employed and the system is simulated via a Monte-Carlo procedure. This 
furnishes the equilibrium properties of the system at its interfaces and 
surfaces as well as in the bulk. 

A. Parametric Investigation 

The following progress has been made during this period: 

1. A utility package has been prepared for the evaluation of the MOLDY 
results. This package consists of several programs and uses the data 
generated by the molecular dynamics program as an input. It has the 

following capabilities: 

(i ) The radial distribution function, virial, density and energy 

profiles across the interface can be obtained. 

(ii) It can perform time series analysis: correlation functions 

(density or velocity auto-correlation) and vibrational frequency spectra 
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of the particles in motion can be calculated. The program has the 
flexibility of calculating components of a vibrational frequency 
spectrum for any number of particles and for any desired subsection of 
the system. It uses a fast Fourier transformation routine and accepts 
velocity components of the particles generated by the molecular 
dynamics. The program requires a relatively small amount of CPU time. 

2. A graphic software package has also been developed which provides visual 
inspection of the whole system as well as of any local region of the 
system. It can be employed for generating a single frame for any 
desired time step or it can be used to produce a movie for the real time 

simulation of individual atomic motions. 

The package also includes a program which can be used to generate 

trajectories of individual particles for quick inspection. 


At the moment a test case is under investigation for the analysis. The 
system consists of two different atomic species to represent an interface 
between a substrate and an adlayer. Figure 1 demonstrates the system 
schematically. It contains 927 particles (with 381 particles (Cu) in the 
adlayer and the rest, 546 particles (Ag), in the substrate). The reduced 

parameters used in the molecular dynamics procedure (representative of a Cu/Ag 
system) were: 
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No periodic boundary condition has been imposed because of the binary nature 
of the svstem. The particles (Ag) within the three bottom layers of the 

substrate were held fixed during the run; but they were permitted to interact 
with the rest of the system contributing to the total potential energy. The 

test system was run up to 1000 time steps. The first 200 steps were used for 

the temperature rescaling process, and the rest employed in the analysis. All 
the thermodynamic quantities were calculated considering only the "core" 
region of the hexagonal slab (to eliminate the surface effect). The core 

included all the atoms within a radius of 2.8 from the center. Figure 

2 demonstrates mobile atoms during the equilibration; and Fig. 3 is a close-up 
of the same region around the interface. Figure 4 shows trajectories of the 
atoms within the top layer of the substrate, and Fig. 5 is a trajectory plot 
for the bottom layer of the absorbatt (which is facing the layer in Fig. 4). 

Figure 6 demonstrates the potential energy profile across the interface 
which is located between the third and the fourth peaks. The difference 
between these peaks Indicates the stronger binding among the adatoms. Figure 
7 is the virial profile similar to Fig. 6. Because of its differential 

character, the virial reflects the mismatch of the sites at the interface. 

The difference in atomic sizes causes the adlayer side to dilate while the 

substrate side contracts. The dilation is exhibited by a negative virial 

while the contraction by an added positive contribution as expected. 

B. Simulation of Systems with Si and 0 Atoms 

For a more quantitative representation of systems involving Si and 0 
atoms we considered a general type of potential energy function which is based 
on two-body and three-body interactions. The most serious difficulty involved 
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in simulation studies of such systems is, in general, due to the potential 
energy function and its ability to represent the system adequately. The 
analytical form of the potential energy function and the evaluation of its 
parameters that are being used are described briefly. 

For a system of N particles the total potential energy ♦ may be 
expressed as: 


N N 


N N N 


i i * i * k 
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( 1 ) 


where, u(r.,r.), u(r. ,r.,r.) ..., u(r. ,r.,...,r^) are the two-three and n- 
ijijK ijn 

body potentials respectively; r. denotes the position of the particle. 

J 

Clearly, the most important term in this expansion is the first term involving 
the two-body potential. Therefore, in the majority of the atomistic 
calculations to date only pairwise additive potentials have been used. In 
this investigation not only two-body but also three-body interactions are 
taken into consideration. Because this expression (Eq. 1) will be used in 
lengthy machine computations, the functional form of 

u(r^,rj) and «^j ) ™st be very simple. The two-body part, therefore, 

is represented by a Mi e-type potential. 
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where ■ Ir^ ’ I • denotes the equilibrium separation and c 

is the energy at r^j » r®. The exponents m and n account for the repulsive 
and attractive terms, respectively. For the intermediate to large separations 
the three-body interactions may be expressed as: 


i 

where, the summation includes all triple multiple interactions resulting from 
the expansion of the third-order interaction energy for three atoms. Each 
term in the sunwnation is expressed as the product of a geometrical factor 
6{r^.rj,r|^) which depends on the relative positions of the three atomic 
nuclei and an interaction constant which depends only on the atomic species 
involved in the interaction. Here, we employ only the triple-dipole 
interaction which has been shown to be the dominant contribution 
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where e., ®., and r.., r.. , r.. represent the angles and the sides of 

1 j k ij 

the triangle formed by the three particles i, j and k. 

Now, combining equations 1 through 5 and neglecting the four-body and 
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higher terms one obtains: 


ORKaNAL P*0* ® 

e<^rtQ OUALITY 


♦ ■ ^2 ♦3 


(6a) 


with 


N N 


1 M 


m 


- m 


1 J 

i j 




(6b) 


and 


1 JL_1_!L (1 3Cose^ Cos6j Cos e^) 


• ; 2:^1: 

® 1 j k 

1 ♦ j k 


* '‘Ik * '‘jk^ 


T 


(6c) 


In this general form the total potential energy ♦ Is a function of the 
parameters, e, r°, m, n and 2j, and the atomic configuration of the 
system. The parameters are assumed to be Independent of the atomic positions 
and the geometry of the system. They depend only on the atomic properties of 
species Involved. “ 

For a binary system of types 1 and 2 [1 = silicon, and 2 s oxygen] with 
the corresponding number of particles Nj and N 2 . Eq. (6a) may be rewritten as: 
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where 
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r.j * 0 (11) 


r,j . 0 (12) 


and 



j < k 




(13) 


r° represents the equilibrium separation between the two species o and 6 
for the two-body part (r°g s *“ 00 ^’ * critical distance parameter 

which Is used as a normalizing factor. * Is Independent of the value of d. 
''aj * 1^0 position of an atom of type o taken 

as the central particle for the summation. The total number of particles In 
the system is given by ■ Nj + N£. 

The parameters (such as ^o0 ’^ae ’®a0 *^a6Y ^o6y^ obtained using 

isolated microclusters and the bulk property data for pure Si and 0, as well 
as data for various SIO 2 crystalline phases. Furthermore, for the evaluation 
of the parameters, a stability criterion was considered that is given as: 

(^ ) . 0 . ,14, 

0 
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For a start, exponents in and n In the two*body term were assumed to be 
equal to 12 and 6, respectively. Values for 

computed directly from the crystallographic data for the corresponding $102 
phases (see Table 1). The parameters silicon were 

calculated from experimental bulk and small cluster data, while for oxygen 
they were obtained from thermodynamic data reported for the O 2 molecules as: 

- 37745 (®K) 
rjj - 2.3516 (A) 

Zjjj - 6.136 X 10^(®K • A®) 

C 22 - 59943 (®K) 
r^2 ‘ 1*2078 (A) 

Z 222 • ~ 0*2623 X 10^(®K • A^) (*) 

Next, we used thermodynamic data (namely cohesive energies) for a-quartz 
and a-cristobal ite. Employing Eqs. (6) and (14), the rest of the potential 
energy parameters were calculated 


*This value is only an approximation at this stage. However, (If needed) It 
may be perfected using bulk SIO 2 data. 
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€j2 • 42531.7 (°K) 

rjg • 1.6259 (A) 

Zii2 • 5.0086 X 10® (“K * A^) 

Zi 22 - 3.9999 x 10® (°K ® A®) 


was 

This set of parameters/ used In Eq. (7) and correctly reproduced the 
cohesive energies for the 6-cr1stoba11te and B-quartz. These experimental 
cohesive energies are tabulated In Table 1. 

For further Investigation of these systems, the following programs have 
been prepared. 

(I) M onte»Car1o Program ; This program Is written to handle systems 
containing up to three different atomic species. It uses Eq. (7) and 
generates appropriate Markov chain configurations of the system In a 
very efficient way. It stores the list of neighbors using the “large 
core memory” (In a COC 7600) that reduces the computation time 
considerably. In partlcjlar, for larger systems. 

(I I ) A Software Package for Structural Characteristics of SIO^ Systems : 

This package contains basically two different programs. The first 

is to calculate the radial distribution function of the systt..-. It 

provides the radial distribution function associated with any portion of 

(R) 

the system; and It Is able to calculate partial (I.e., 

8(SI-0) '>'■ 3(0-0)> “**’ g(R) of any SIO 2 system. The second 

program Is to calculate angular distributions for the whole system or 
for any subsection of the system. In particular. It Is designed to 
calculate the distributions for angles such as (S1-0-S1) and (0-S1-0) 
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between ihe neighboring atoirs. 

Initial simulations were designed. In particular, to check the 
adequ^%.y of the computational methods and the potential energy 
function. First, the bulk and the surface structures of crystalline SI 
and two SIO 2 forn»^ (namely, a-cr1stobal1te and o-quart^l were taken 
Into consideration. In all cases the very same set of potential 
parameters were employed. 

The results were found to be very encouraging. They Indicate that 
(1) the type of the potential energy function Is quite appropriate for 
these systems, and (11) despite relatively long three>body Interaction 
calculations, the simulation procedure (MonU-Carlo) Is able to produce 
the equilibration (for systems of reasonable sizes) within manageable 
computation times. 

1. Bulk Silic on 

The crystalline structure of silicon was generated in a diamond cubic 

3 

structure for T ■ 300®K corresponding to an atomic volume of ?0.02 A (exper1> 
mental data were taken from reference 8). The system contained 216 SI 
atoms. To mimic an infinite bulk a three-dimensional periodic boundary 
condition (PBC) was employed in accordance with the experimental density. The 
result of an equilibrated Monte-Carlo run indicated that the crystalline 
structure of the system is well maintained and the stability criterion was 
observed. 

For a more stringent test one of the high pressure forms of Si (in the 
white tin structure) was considered. This allotropic form is much denser than 
the regular diamond cubic structure of Si. For both for*ms we investigated the 
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variation of the potential energy as a function of volume at the static 
limit. The result is shown in Fig. 8. Curves intercept each other at 
V • 14.95 which indicates that, for V < 14.95, the white-tin structure 
would be energetically more favorable than the diamond cubic structure. This 
outcome agrees quite well with an atomic volume of 14.95 which was 
measured experimentally for the white-tin structure of Si [8]. Unfortunately, 
in this experimental report the pressure to which this volume corresponds was 
not stated. 

2. Si (100) Su rface: 

In this case, we first generated an ideal (1 x 1) surface which was 
obtained simply from truncating the bulk crystal. The same lattice dimensions 
were used as above. This surface structure was taken as an initial 
configuration. The top layer contained 32 Si atoms out of 215 which 

constituted the system. For the simulation, PBC was employed in two 
directions (i.e., in X- and Y-directions) leaving the free surface in the Z- 
direction. In the relaxation process only the top layer Si atoms were 

considered. The rest of the atoms were held fixed in their original lattice 
sites but they contributed fully to the interaction energies. Figure 9 
displays the surface configuration after the equilibration.* The relaxation 
produced a (2x2) structure consistent with the recently proposed models 
[9,10]. 


*In all '•ases, iterations were carried out beyond the complete equilibration 
point wh' :.i was monitored by the variation of the total energy. 
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Bulk g-Cristobalite and a-Quartz 


! Jk a-cristobalite and o-quartz were simulated to 50"K, 300*K and 
2500®K. PBC was imposed in three dim6n«5ions in accordance with the 
experimental densities at the corresponding temperatures and species [2]. At 
T » 50®K and 300®K, the equilibrated systems maintained their crystalline 
structures quite well. The configuration of these "relaxed" structures was 
analyzed using the radial distribution function (RDF) and the angular 
distf'bution function (ADF) calculated from the final atomic positions. For 
both species the RDF produced peaks in the correct positions compared to the 
ideal crystals. Figure 10 displays the RDF for o-cristobal ite as an example 
at T = 300®K for SiO, 0-0 and Si-Si cases as well as the total RDF. Figure 11 
shows the corresponding RDF's for an ideal case. Calculated ADF's for both 
a-cristobalite and o-quartz produced Si -0-Si and O-Si-0 angles consistent 
with experimental results [2]. The distribution of O-Si-0 angles averaged 
around 109* for both cases confirming the inherent tetrahedral structure of 
the system. The spread In this distribution at T « 300*K was approximately + 
10*. The distribution of Si -0-Si angles at 300*K for a-cristobalite varies 
from 139* to 169* with an average of 154*, and, for a-quartz, it varies from 
132° to 165* with an average of 147*. Figure 12 displays the ADF's for 
a-cristobalite at T = 300*K and the corresponding ideal crystal case (i.e., 
before equilibration) is shown in Fig. 1?. 

To test the high temperature behaviors of the system a simulation run was 
carried out at T = 2595*K. The PBC imposed on the system was chosen such that 
it produced the correct density (1.929 gr/cc) at this temperature. RDF 
obtained after the complete equilibration produced more diffused peaks than 
the lower temperature cases which is indicative of a molten state (see Fig. 
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14). However, the general feature of the total RDF was still preserved. For 
the ADF case also we obtained peaks somewhat more spread. The distribution 
for the Si -0-Si angles varies from 145® to 180® averaging around 165 
degrees. This indicates a shift to larger angle region with increased 

temperature, as expected. The distribution of O-Si-0 angles, on the other 
hand, was found well around 109® region (with a spread of + 15®). This 
confirms that the stability of the tetrahedral units, even at this high 
temperature, is still maintained (see Fig. 15). 

These results, in particular those for 50®K and 300®K, are very important 
and indicate that the present semi -empirical potential function is quite 
adequate and it is able to reproduce stable crystalline structure of 
a-cristobal ite and a-quartz satisfactorily. Of course, the most 

interesting part is the preservation of the crystalline structures in the 

Monte-Carlo scheme at finite temperatures. At the present, there is no simple 

potential function in the literature that can be used in a similar computer 

study for crystalli.ie Si02 forms. The only model potential function which can 

be found in the literature is for the amorphous Si02« It is based on an ionic 
model and employs a two-body Born-Mayer-Huggins-type potential function [3- 
6]. This potential, however, fails to reproduce an acceptable crystalline 
SiO^ form as a stable state. 

These present results are particularly important because they clearly 
demonstrate the role played by the three-body forces in the stability of 
crystalline Si 02 forms. 
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4. Surfaces of a-Cristoba11te and a-Quartz 

Surface structures for the (001) plane of a-cristobalite and (OOCi) 
plane of a-quartz were taken Into consideration. These surfaces were first 
generated by truncating the corresponding ideal crystals as in the case of 
Si (100) explained above. Both of these surfaces are 0-rich surfaces. Two 0 
atoms from each surface tetrahedron are exposed. Figure 16 demonstates top 
and side views of the surfaces in their initial configurations. Both systems 
contained 216 atoms. For a-cristobalite, only the top 54 atoms and, for 
a-quartz, only the top 72 atoms were included in the relaxation procedure. 
The rest of the atoms were held fixed in their lattice positions, but they 
contributed fully in the energy calculations. PBC was applied in two 
dimensions to provide a continuous surface layer. Satisfactory equilibration 
was reached after approximately 1500 iterations. Top and side views of the 
equilibrated surfaces are shown in Fig. 17. Both surfaces exhibited similar 
patterns. The top layer 0 atoms of the neighboring tetrahedra were paired and 
formed so called "peroxide bridges," [11] with an average 0-0 distance of 
1.2 A. Figure 18 displays the RDF for the relaxed a-cristobalite surface 

region. The first peak is due to the peroxide bridge which does not exist in 
the RDF's of bulk structures. All other peaks are approximately in their bulk 
positions. On these reconstructed surfaces one can distinguish three 
different types of tetrahedral units (i ) with no peroxide-bridge, (ii) with 
only one peroxide-bndge and (iii) with two peroxide bridges connected to both 
sides. At the moment the energetics and the exact geometry of the different 
tetrahedral unit' are not determined. 

This preliminary test clearly indicates that, in the study of surfaces 
the role played vy the three-body interactions is very important, and the 
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present potential energy function is well suited for such investigations. 
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TABLE 1 


g-Crist. a-quartz 6-Crm. B-quartz 


R(Si-Si) 

3.0767 

3. 

R(Si-O) 

1.60 

1. 

R(O-O) 

2.60 

2. 

a"* 

4.0575 

4. 

a*" 

hz 

9722.33 

8328. 

a*" 

®22 

41.3494 

36. 

a*" 

®21 

4861.17 

4164. 


5.2401 

5. 

®12 

202.442 

188. 

a" 

®22 

19.5795 

20. 

a" 

®21 

101.221 

94. 

^U1 

5.7840 

8. 

r 

112 

-191.313 

-117. 

^122 

174.512 

195. 

T222 

62.0504 

76. 

''’221 

174.555 

195. 

"^211 

-47.8405 

-29. 

Cohesive 

Energy 

(kcal/mole) 

-412.14 

-412. 


0579 

3.0880 

3.0982 

60 

1.55 

1.63 

63 

2.53 

2.60 

1148 

4.0388 

4.1048 

37 

9796.98 

9870.78 

7629 

42.0266 

48.1561 

88 

4898.46 

4929.74 

6759 

5.0999 

5.6325 

962 

202.892 

204.702 

2410 

19.3041 

21.5124 

4946 

101.444 

102.302 

3242 

4.94088 

8.0688 

663 

-202.276 

-175.685 

381 

167.158 

195.231 

5891 

57.5318 

73.9876 

427 

167.131 

195.237 

4278 

-50.5788 

-43.9041 

540 

-417.43 

-417.85 
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FIGURE CAPTIONS 


1. Schematic representation of atomic system for interaction between copper 
and silver. 

2. Simulation of Cu movement on Ag(lll). 

3. Close-up view of Fig. 2. 

4. Trajectories of the atoms within the top layer of the Ag{lll) substrate 

at the Cu/Ag interface. 

5. Trajectories of the atoms within the bottom layer of the Cu adsorbate at 
the Cu/Ag interface. 

6. Potential energy profile across the interface illustrating the stronger 
binding among the adatoms. 

7. Vi rial profile across the interface reflecting the mismatch of the sites 
at the interface (adlayer side dilates while the substrate side 
contracts) . 

8. Cohesive energy of silicon as a function of its ''olume. Curve 1 is for 
diamond cubic structure and curve 2 is for white tin strucutre. Energy 
values are reduced by the two-body energy parameter, e, and the volume 
is in A^. 

9. Reconstructed ‘‘i(lOO) surface. Largest symbols represent top layer Si 
atoms. Smaller symbols are 2nd and 3rd layer atoms. Only top layer Si 
atoms were relaxed. 

10. Radial distribution function for bulk o-cristobal ite at T « SOO^K, 
after the equilibration, (a), (b), (c) and (d) display distribution of 
Si-0, 0-0, Si-Si bonds and the total RDF, respectively. 
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11. Radial distribution function for bulk a-cr1stoba11te, Initial state, 
(a), (b), (c) and (d) are the distribution of SI-0, 0-0, SI-SI bonds and 
the total RDF, respectively. 

12. Percentage distribution of bond angles for a-cr1stobal1te at 300**K 
after equilibration, (a) and (b) exhibit distributions for 0-S1-0 and 
SI -0-SI angles respectively. 

13. Percentage distribution of bond angles for o-cristobalite at the 

initial state, (a) and (b) display distributions for O-Si-0 and Si-O-Si 
angles, respectively. 

14. The radial distribution function at T * 2595*K. (a), (b), (c) and (d) 

display distributions for Si-0, 0-0, Si-Si bonds and the total ROF, 
respectively. 

15. The distribution of bond angles at T » 2595“K. (a) and (b) are the 

distribution of O-Si-0 and Si-O-Si angles, respectively. 

16. a-cristobalite (001) surface structure at the initial state, (a) Top 

view, (b) and (c) are the side views. a-quartz (0001) surface 
structure at the initial state, (d) top view, (e) and (f) display side 
views. 

17. Equilibrated structure for o-cristobal ite at 300®K. (a) Top view (b) 

and (c) are side views; and equilibrated structure for a-quartz at 

SOO^K (d) top view (e) and (f) are the side views. 

18. The total radial distribution function for a-cristobalite at 300°K 
after equilibration. Only surface region atoms are included. 
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FIGURE 0 - 3 
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In many areas related to nucleation, surface reconstruction, surface 
transport, adsorption and desorption, the detailed mechanism associated with 
the migration of clusters Is highly desired. Several models based on 
experiments of dimer and trimers motions have been reported [1-3], Detailed 
quantitative data for cluster motion on substrate surfaces are limited to 
Vl(211) and W(llO) surfaces. On these surfaces much of the work on clusters 
has been concerned with their structural stability [3,4]. Trends In the 
dependence of cluster structure as well as dimer binding-energy on chemical 
Identity are emerging, but the data are still limited. No computer simulation 
study has been yet reported for an atomic level analysis of the mechanism 
associated with the motion of clusters at substrate surfaces. 

In this part of the study the diffusion process for dimers (such as Pt 2 
and AU 2 ) on a Pt(llO) substrate surface was simulated employing a molecular 
dynamics technique based on Lennard-Jones type two-body Interactions. Because 
of the stability problems caused by this potential energy function, here, we 
considered a Pt(llO) surface as a substrate whose atomic arrangement Is 
similar to a W(2H) surface. Calculations were carried out at Tp^ « 0.3; 
where, Tp^ denotes the reduced temperature of Pt. For gold this reduced 
temperature corresponds to .46. The reduced time step was taken as 

dt* < 0.01. (in real units the time step Is on the order of 10"^^ second). 
The reduced lattice constant, a*, of the substrate was calculated in two 
ways. (i) Using the density and the thermal expansivity of Pt. These were 
taken from the literature as d 2 o » 21.45 gm/cc and * 8.9 x 10"®/deg [5], 
respectively. The melting point of Pt (2045*K) was also used to scale the 
temperature of the system to give ^JJielting * which is the case for 
Lennard-Jones systems. Considering the Lennard-Jones parameters for Pt(cp^ 
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« 7910«K and Op^ « 2.542 A [6]) we calculated the a* as 1.5518, at T* » 
0.3. (i1) Using the atomic volumes of rare gas solids, at t* ■ 0.3 the 

reported values for the molar volumes of Ar and Xe are 22.93 and 35.834 
cm^/mole, respectively [7]. Accordingly, the reduced lattice constants were 
obtained for Ar as 1.5706 and for Xe as 1.5570. These were calculated 

employing the Lennard-Jones parameters reported in reference [8]. The value 
of a* obtained from the Pt data produced a somewhat compressed substrate 
structure, and the i* obtained from the Ar data produced a dilation that was 
monitored by negative and positive virial values during the molecular dynamics 
run. In the study of a Pt/Pt(110) system even a slight compression of the 
substrate (i.e., a smaller a*) was found to have considerable effect on the 
exchange process. Compression encourages the cross-channel diffusion while a 
dilation favors channel diffusion. In the rest of this study we employed a* • 
1.559 which is an averaged value for the reduced lattice constant calculated 
from Pt, Ar and Xe data. 

To simulate the substrate, a rectangular slab of atoms was generated with 
periodic boundary conditions ir the x and y directions with the (110) surface 
normal to the - direction. The dimers (Au 2 and Pt 2 ) were positioned on the 
(110) surface. For both cases, the atoms were placed on the adjacent channels 
to study the dimer motion. For the dynamic calculations only the upper 192 
substrate atoms and two adatoms were mobile. The rest of the system; i.e., 
the lower 192 substrate atoms were kept inmobile at their lattice positions. 
However, these fixed atoms were taken into account in the total energy and 
force calculations. Velocities of the mobile atoms were initialized with a 
Maxwellian distribution to the temperature of the run. For equilibration, 
velocities were rescaled every tenth step for 500 steps. An additional 3000 
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or more steps were used to observe diffusive behavior of the surface atoms. 
Figure 1 displays the trajectories for the Au dimer and the top layer of the 
Pt substrate atoms for 3400 steps. Circles are centered on the final atomic 
positions. These trajectories indicate that the dimer has been associated on 
the Pt(llO) surface. One of the Au atoms moved two lattice sites apart while 
the other Au remained localized. In a different run for the Au 2 /Pt case, one 
of the Au atoms formed a split interstitial and a temporary exchange which is 
similar to the single Au atom diffusion case [9]. 

Figure 2 displays the trajectories for the Pt dimer along with the top Pt 
atoms of the substrate. In this case, the motion of the dimer remained quite 
localized up to 10,000 time steps. Only 0«state or 1-state configurations 
were observed [1-2]. In another run for Pt2/Pt(110), considering a slightly 
smaller a* (1.5518), it was found that one of the Pt adatoms experienced an 
exchange with the wall atom (cross-channel diffusion). From this 
investigation it appears that the exchange process is quite sensitive to the 
value of a*. This is an important outcome and has never been investigated 
before. 

Two important factors may contribute to the dissociation of the Au 
dimer, (i) The reduced temperature, T*, of Au is higher than the T* of Pt, 
and (ii) the ?t-Pt interaction is considerably stronger than the the Au-Au 
interaction. In this study, no many-body effects were taken into account. 
However, in the dimer dissociation process, many-body interactions may play an 
important role. The effect of three-body interactions along with the effect 
of T* on the dissociation process of dimers on substrates remain to be 
investigated. 
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FIGURE CAPTIONS 


1. Trajectories for a Au dimer and the top layer of Pt substrate atoms. 

2. Trajectories for a Pt dimer and the top layer of Pt substrate atoms. 
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CRACK FORMATION AND PROPAGATION PHENOMENA 



A. 


Introduction: 


The mechanical behavior of materials can be specified by macroscopic 
theories on the basis of a fcM material constants which provide an accurate 
description for the responses of materials to forces. However, these theories 
don't provide a microscopic level of understanding of these basic mechanisms, 
such as those underlying plastic deformation. In general, the theories for 
strength of materials, elasticity and plasticity lose much of their power when 
the structure of a material becomes an important consideration and the 
material can no longer be considered a homogeneous medium [1]. 

It is well recognized that the relationship between mechanical behavior 
and microscopic structure of materials is very important. When mechanical 
behavior is understood in terms of microscopic processes, it is often possible 
to improve the mechanical properties of a material. A microscopic description 
of materials, as opposed to macroscopic theories, cannot be achieved by 
properly defining a few material constants. Instead, the system must be 
described in terms of interatomic forces and the coordinates of the particles 
which constitute the material. 

Several studies based on the computer simulation of the mechanical 
behavior of materials exist. These studies have provided a better under- 
standing of various mechanisms [2,3], sue* as diffusion, crack propagation, 
dislocation motion, plasti*" flow, etc-, at the atomistic level. 

Material failure is generally caused by fracture which follows yielding 
or plastic deformation. Slip is the simplest and most common example of 
plas.ic flow encountered prior to the ductile frature of materials. In this 
investigation, the effect of a uniaxial load exerted on a two-dimensional 
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microscopic slab was analyzed by a computer simulation technique. Even though 
the discrete nature of the model has limited the size of the sample, the study 
has produced illuminating results which agree well with the macroscopic 
theories and provide an atomistic level of understanding of slip and the early 
stages of crack formation. 


B. Computational Model and Procedure 

The irodel consisted of a two-dimensional triangular lattice which 
contained 400 or more particles. Each particle was treated discretely. The 
particles were assumed to be interacting via a Lennard-Jones type two-body 
potential function which is given by: 


U(r,j) 



( 1 ) 


where is the pair potential energy between particles i and j separated 

by r^j. The equilibrium distance is denoted by Cq and e is the interaction 

energy at r^^ = r^. Many-body interactions were neglected in this 

investigation. All particle neighbors up to 3.5 r^ were included in the 

energy force calculations. This procedure insures that all neighbors up to 
the fifth-nearest neighbor will be included in the calculations. In the 
calculation procedure, periodic boundary conditions were not employed. This 
allowed us to examine the surface reconstruction and the formation of edge 
dislocations during the elongation process. 

Initially, the system was generated in a rectangular shape in its 
equilibrium configuration. Then, the system was elongated in a stepwise 

fashion by imposing a uniaxial load via small increments. After each 
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Incremental elongation, the system was allowed to fully equu'trate during a 
relaxation period. All of the mobile particles (defined below) were fully 
equilibrated after every incremental elongation. A force minimization 
technique was employed to do this. The force acting on each particle was 
calculated; then, the particle was moved in the force direction until the 
resulting force acting on it became virtually zero. This procedure was 
repeated sequentially for every moving particle up to the complete 
equilibration state. 

The crystal contained two types of particle: (i) Mobile particles: the 
majority of the particles were treated as mobile particle. These particles 
were located in the central portion of the lattice as illustrated in Fig. 1. 
During the elongation procedure, their positions were displaced uniformly in 
the load direction according to a prescribed incremental elongation. All 
mobile particles were fully equilibrated during the subsequent relaxation 
procedure, (ii) Fixed particles: these particles constituted the rigid edges 
and were held immobile during the relaxation period, i.e., the relative 
positions of the fixed particles with respect to each other remained unchanged 
within each particular rigid edge. However, these particles were permitted to 
interact with the mobile particles, and therefore, contributed to the total 
energy as well as to the force calculations. 

C. R esults 

1 . Two Dimensional Perfect Triangular Latti ce 

A perfect triangular lattice (i.e., the basal plane of an hep crystal 
[4] under a tensile load was studied. The tensile load was applied to three 
different orientations of the lattice. In the first case, the load was 
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applied in the [1T2] directions; in the second case in the [OlT] direction; 
and in the third case, the load was imposed in the [F14] direction (see Fig. 
2). In all three cases, the two dimensional slab, because of the rigid edges, 
experienced a tensile deformation with a "constraint” [1]. 

Small incremental elongations* were Imposed on the system which resulted 
in a progressive increase in the total energy until the first slip event 
occurred when, as shown on Figs. 3(a)-(c), the energy decreased 
dramatically. Continued elongation of the lattice resulted In additional 
occurrences of slip. Each slip was accompanied by a rapid decrease in the 
total energy as shown in Fig. 3(a), for the [112] case. Up to the 
occurrence of the first slip event, the symmetry of the lattice was well 
maintained, i.e., only a small amount of lateral contraction was detected at 
larger elongations. However, no attempt was made to estimate the 
corresponding Poisson's ratios. 

In each case, slip always occurred in the close-packed directions as 
predicted by macroscopic theories. Since close-packed directions were 
differently oriented with respect to the tensile axis, slip Initiated at 
different elongations. From the variation of the calculated total potential 
energy (as a function of elongation), the "average" uniaxial stress and the 
resolved shear stresses were calculated as described below, and tabulated in 
Table 1. 


The incremental elongation which was used (around the slip formation) 

corresponded to an imposed fractional strain of 0.001. The magnitude of the 

incremental elongation had some effect on the results. A test study indicated 

that results were virtually unchanged for incremental elongations up to ne * 
0.1. However, larger Ae values caused the system to undergo a sudden 

brittle-type fracture. 
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from; 


The resolved shear stress, along the slip plane was estimated [1] 


F 

m (2) 

A 

where F/A denotes the average uniaxial stress; F is the uniaxial load, and A 
is the cross-sectional area The variable m is the Schmid factor and is given 
by: 


m = cos ♦ • cos X (3) 

where, ♦ is the angle between the "normal" to the slip plane and the tensile 
axis and x is the angle which the slip direction makes with the tensile axis 
as illustrated in Fig. 4. 

The value of F was calculated by numerically differentiating the total 
potential energy, E, of the system with respect to the elongation: 

dE 

F = - (4) 

dL 

where L represents the length. As an example. Fig. 5 demonstrates for the 
[Tl2] case, the stress-strain curve, which is the variation of F/A as a 
function of (L-Lq)/Lq, where Lq is the initial lenath. In Table 1, the values 
of F/A were calculated using the values of E and L just before the occurrence 
of slip. In the present model, A is simply the width of the crystal. 
Throughout this investigation, the units of F and A were normalized with 
respect to c and r^ of the Lennard-Oones potential, (i.e. F is given in 


units of E/r^). Table 1 also contains the values for the tensile stresses, 
0 ^, normal to the slip plane which were calculated as: 

^ 2 

o_ » _ • (cos • (5) 

' A 


The values of Tj^, calculated from Eq. (2) and listed in Table 1 for 
three different load directions are relatively close to each other. 
Obviously, the magnitude of the stress acting normal to the slip plane has an 
effect on the process of slip in the present two-dimensional model. For 
example, for elongations in the close-packed [112] direction, the relatively 
low resolved shear stress is associated with a higher tensile stress (normal 
to the slip plane) which tends to separate the neighboring close-packed rows 
facilitating the slip. 

Formation of slip events associated with the elongation process are 
demonstrated in Figs. 6(a) to 6(e) for strain applied in the close-packed 
[TT2] direction, in Figs. 7(a) to 7(e) for the perpendicular [OlT] 
direction and in Figs 8(a) to 8(b) for the [F41] direction. In all three 
cases, slip occurred along the most densely packed rows as anticipated from 
the macroscDpic theories. The triangular lattice considered here has three 
close-packed rows indicated by axes a^ , ^3 2. 

For strain applied in the [TT2] direction (i.e., along the a 3 axis) the 
other two close-packed rows (i.e., a^ and a 2 ) form a 60® angle with the 
tensile axis. Figures 6(a) and 6(b) represent the fully equilibrated system 
with an average linear strain, e, of 0.0455 and 0.0974 respectively. Figures 
6(c) and d demonstrate the initiation and the completion of the first slip 
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(which is in the CT2T] direction) after an incremental elongation corres- 
ponding to e * 0.0989. Figure 6(c) is not fully equilibrated energy wise, but 
it shows the initiation of the slip via an edge dislocation type process. 
Figure 6(e) represents the system with e ■ 0.1499 after several slip events 
have occurred. Because of the rigid edges, the system was not permitted to 
deform freely by uniform glide. Therefore, the slip planes rotated toward the 
tensile axis. This behavior is consistent with experiments [1]. Results of 
similar calculations indicated that slip along either aj or a 2 occur with the 
same probability. 

For the load applied in the [OlT] direction, one of the close-packed 
rows (i.e., the a^ axis) is perpendicular to the load direction. In this 
case, the other two close-packed rows, namely, a 2 and a 3 both form 30® angles 
with the tensile axis. Figure 7(a) shows this system under a tensile load in 
the [OlT] direction with e ■ 0.1270. This was the last elongation-relaxation 
step for the system before slip occurred. The next incremental elongation (e 
= 0.1291) caused several concurrent slips to occur and the total energy of the 
system was lowered about 0.25 e. Figure 7(e) shows the system after it has 

fully relaxed. Figures 7(b), 7(c) and 7(d) represent intermediate states 
during the relaxation procedure. They demonstrate clearly the formation and 
the motion of edge dislocations associated with slip and vacancy formation 
[1,2]. Figure 7(d) shows a dipole where dislocations have arrested each 

other. Further elongation caused the system to break in two parts where the 
vacancy was formed (indicated by an arrow in Fig. 7(e). 

In the case of the tensile load applied in the [541] direction, the 
three close-packed rows, a^, a 2 * and a^, make 19®, 41® and 79® angles with 
the tensile axis, respectively. Figure 8(a) demonstrates the system with an 
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average linear strain of 0.0786 before the slip. Further elongation 
corresponding to e * 0.810 caused the slip to occur along the ag *xis that 
is shown in Fig. 8(b). The results indicate that the probability of 
occurrence of a slip event increases as the angle formed between the close- 
packed row and the tensile axis approaches 45°. This is consistent with 
experimental findings [1]. 

The results of this study indicate that perfect lattice structures 
exhibit relatively high resistance to Jl'p and deformation. Because of the 
two-dimensional character of the model and the neglect of the kinetic energy 
(i.e., temperature effects), the present results should only be used for 
qualitati/e comparisons with experimental data. However, the absolute values 
of the theoretical results necessary to produce slip are in reasonable agree- 
ment with the yield strains obtained from experimental tests on whiskers [6], 

In the tensile strength experiments (which are believed to be nearly 
perfect crystals) a large variety of materials are found to have quite high 
yield strains [5,6]; in fact, comparable to those obtained in this study. The 
strongest whiskers are the smallest in size. As the diameter and length is 
increased, the strength of the whiskers decreases considerably. It is 
postulated that this decrease in strength is due to defects which are formed 
accidentally during the growth of the whiskers [5]. In the following section 
the influence of points defects will be investigated. 

2. Two-Dimensional Lattice with Defects 

The effect of placing a point defect near one of the surfaces in the 
two-dimensional model is shown in Fig. 9, As illustrated, the tensile load is 
applied in the (TT2] direction. Figure 9(a) shows the model at equilibrium 


75 



before the load Is applied (e ■ 0.00). The point defect remains in its 
original position until the strain reaches approximately 0.066 when, as shown 
in Fig. 9(b), the defect moves from the interior of the lattice to the surface 
as expected [7]. The vacancy remains fixed on the surface as the strain is 
increased until slip occurs at e ■ 0.077 (Fig. 9(c)). The slip in the [T?T1 
direcUon was initiated by the vacancy. Additional simulations in which the 
defect position was varied while maintaining the distance between *'he defect 
and the surface show results similar to those shown in Fig. 9. In each of 
these tests, the defect remains in its original position until the strain is 
approximately 0.066 when the defect moves to the surface creatin a vacancy 
there. The vacancy remains on the surface until slip occurs. The onset of 
slip always occurred at e ■ 0.077, however, in some of the tests, the slip 
occurred in the [2TT] direction. This is to be expected since slip along 
both directions are equally probable. The strain (e » 0.77) at yhich slip 
initially occurred is less than the strain of 0.098 at which the "perfect" 
crystal showed sli,;. This reduction in strain, necessary to produce slip, is 
indicative of the role that defects have in weakening materials. 

Figures 10(a) - 10(c) show the results of another simulation in which the 
point defect is placed near one of the surfaces. The tensile load in this 
test is applied in the [OlT] direction. The results are similar to those 
described above. In this test, the defect moves to the surface when the 
strain reaches approximately 0.047 (see Figs. 10(a) and 10(b)). This is 
expected because in this configuration, the particle separation which prevents 
the defect from being filled is parallel to the tensile direction. Slip, 
however, does not occur until the strain reaches 0.1107 (Fig. 10(c)). As 
previously discussed, these results were very repeatable with slip occurring 
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along both of the equally probable close-packed planes (i.e., [TTZ] and 
CT2T3 directions). 

It is not clear from observing Fig. 10(c) that the vacancy initiated the 
slip. However. Figs. ll(a)-lUd), which represent intermediate states during 
the relaxation procedure, show how the vacancy is related to the occurrence of 
slip. The slip which occurs in the [T2T] direction is triggered by the 
vacancy and progresses via a dislocation motion which is clearly visiole in 
Figs. ll(a)-(d). Near completion of the slip process (Fig. 11(d)) however, 
the slip direction changes to CTT2]. The completed slip process is shown in 
Fig. 10(c). 

Figures 12 (a) -(c) show the results of a test to determine the effect of 
placing the point defect deeper in the lattice. Figure 12(a) represents the 
model with e > 0.0. In contrast to the results in which the defect is placed 
near the surface, the vacancy does not move to the surface as the load is 
applied. Instead, the defect remains at its original position until slip 
occurs (Fig. 12(c)). In the tests performed, slip always occurred along a 
plane adjacent to that containing the void. An intermediate state (Fig. 
12(b)) shows that the defect moves as the slip process takes place. Slip was 
detected at a strain of 0.077 which is identical to that observed when the 
defect was initially placed closer to the surface. 

0. Conclusions : 

The discrete model used in this investigation produced results that are 
consistent with those from macroscopic theory [1]. The results described 
herein indicate that slip events, which are known to be the simplest type of 
plastic deformation in crystalline bodies, occur predominantly on rows with 
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higher density atoms (i.e. along the close-packed rows). In addition the 
calculations for perfect 2-d1mens1ona1 triangular crystals with uniaxial loads 
Imposed In the [OOT], [112], and [514] directions showed that the rule of 
maximum resolved shear was observed [1]. These simulation results also 

clearly illustrate the involvement of dislocations in the slip formation 
process as anticipated from macroscopic theories. 

The results of tests to determine the impact of varying the incremental 
elongation ratio showed that the general characteristics of ductile-type slip 
formation are not affected by the magnitude of small incremental elongations 
corresponding to fractional strains smaller than 0.01. However, larger 

increments In elongation cause the material to undergo a sudden brittle-type 
fracture. These tests also showed that, as the single crystal deforms under 
the uniaxial tension, the slip planes tend to rotate toward the tensile 
axis. This is due to the rigid edges which force the grips to remain in line. 

In all cases studied, the systems with point defects experienced slip 
formation at smaller strains than the corresponding perfect crystals. 
Vacancies located near the surface move to the surface before the slip 

occurred. However, vacancies In the Interior regions move to the surface only 
during the slip process. 

The present Investigation based on energy minimization has limited 
applicability. It Is used here chiefly to test the model, however, the 

results were quite Illuminating In providing atomistic level understanding of 
various phenomena related to plastic flow. Furthermore, caution should also 
be used in Interpreting these results due to the two-dimensional character uf 
the model. Further studies are In progress, in these laboratories, on 30 
systems and employing a molecular dynamics technique to analyze the 
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analyze the temperature and time dependent behaviors of the various processes 
which are closely related to slip as well as crack Initiation phenomena. 
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TABLE 1 


Uniaxial Load Direction 



[II2] 

lOlIl 

(541) 

Angle between the slip plane and the load 
direction, X (in degrees) 

60 

30 

41 

Total number of particles in the system 

400 

414 

420 

Number of moving particles 

320 

324 

340 

Incremental elongation imposed during the 
slip formation (percent) 

0.1% 

0.1% 

C.1% 

Initial cross section. A, (in units of r ) 
Because of the two-dimensional nature of ° 
the model A represents the width 

7.8 

8.5 

8.9 

Total linear strain at the first slip event 
(dimensionless) 

0.C98 

0.128 

0.080 

Uniaxial averaged stress before the first 
slip, F/A (in units of t/r^) 

0.354 

0.457 

0.336 

Critical resolved shear stress, t. 

2 ^ 
(in units of e/r^ ) 

0.153 

0.198 

0.166 

Tensile stress normal to the slip plane, o- 

0.265 

0.115 

0.145 


2 

(in units of ) 
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FIGURE CAPTIONS 


1. Schematic representation of the two-dimensional model 

2. Direction indices in a basal plane of an hexagonal crystal (see Ref. 4). 

3. Variation in total potential energy as a 'unction of the elongation. 

The length is g'ven by L/Lq, where L is the length at any elongation and 
Lq is the initial length. Tensile load applied in (a) CTT2] 
direction, (b) [OIT] direction, and (c) [F41] direction. 

4. Illustration of the geometry of slip. F and represent the tensile 
force and the resolved shear stress, respectively. ♦ is the angle 
between the "normal" to the slip plane and the tensile axis, and \ is 
the angle which the slip direction makes with the tensile axis. 

5. Stress-strain curve for the [1T2] case. The points were calculated by 
numerical differentiation of the energy E with respect to L-Lo)/Lq* 
(Stress is given in units of 

6. Perfect triangular lattice under a uniaxial load in the [TT2] 

direction. 

(a) strain = 0.0455; system fully equilibrated. 

(b) strain = 0.0974; system fully equilibrated before the slip. 

(c) strain = 0.0989; system partially equilibrated. Initiation of the 
slip is clearly visible. 

(d) sirain = 0.989; system fully equilibrated. 

(e) strain = 0.1499; system fully equilibrated. This shows the lattice 
after sever? 1 slip events have occurred. 

7. Perfect triangular lattice under a uniaxial load in the [OlT] 

di rection. 

(a) strain = 0.1270; system fully equilibrated just prior to the 
occurrence of slip. 


83 


(b), (c), (d) - strain « 0.1291; these partial degrees of equilibration 
show the progression of the formation of slip via dislocation type 
motions for three different sequences of the relaxation process. 

(e) strain = 0.1291; system fully equilibrated. Further elongation 
causes the lattice to break. The lattice failed where the vacancy, 
indicated by an arrow, has formed. 

8. Perfect triangular lattice under a uniaxial load in [F41] direction. 

(a) strain * 0.0786; system fully equilibrated prior to slip. 

(b) strain » 0.0810; system fully equilibrated after the first slip 
event occurred. 

9. Triangular lattice with defect (near the surface) under a uniaxial load 
in the [TT2] direction. 

(a) strain = 0.00; system fully equilibrated. 

(b) strain = 0.0644; system fully equilibrated. 

(c) strain = C.0770: system fully equilibrated. 

10. Triangular lattice with defect (near the surface) under a uniaxial load 
in the [OlTl direction. 

(a) strain = 0.00; system fully equilibrated. 

(b) strain = 0.0473; system fully equilibrated. 

(c) strain = 0.1107; system fully equilibrated. 

11. Intermediate steps in the equilibration process of a triangular lattice 
with a defect (near the surface) under a uniaxial load in the [OlT] 
direction. The equilibration process required approximately 300 
iteration steps to reach equilibrium (Fig. 10(c)). 

(a) Step 163 

(b) Step 180 

(c) St . 0 198 

(d) Step 218 
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12. triangular lattice with defect under a uniaxial load in the [TT2] 
direction, 

(a) strain * 0.00; system fully equilibrated. 

(b) strain * 0.0770; system partially equilibrated. 

(c) strain « 0.0770; system fully equilibrated. 
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